Introduction

Our philosophy in this iteration of the analysis is that we can simply treat the number of transects as a covariate to be controlled statistically. This avoids the need to bootstrap and throw away a lot of data, but we pay a price in terms of model identifiability due to the collinearity between elevation and transect number.

Prepare environment

Load data

Calculate network/group-level metrics

This is straightforward because none of these metrics involve independent abundance data (i.e. our survey data). We will focus on the network-level metrics of connectance and NODF (nestedness) and the group-level metrics of niche overlap and C score (both are measures of dietary overlap / partitioning)

Temperature ~ elevation

The question has been raised of whether we might learn more by using temperature instead of elevation as a predictor variable. We see here, though, that total degree-day accumulation and elevation are almost perfectly correlated, so I think there is little to be gained by using temperature as an explanatory variable. We do, however, see that 2011 and 2012 were substantially warmer than 2010.

Extract species-level metrics

How well does mean proportional marginal total predict mean proportional flower cover?

It’s actually not so bad.

Since we have curvy relationships and nested variables, we will analyze our data using hierarchical generalized additive models (HGAMs)

I think there are two ways we could approach this.

  1. If our goal is simply to explain network metrics, we can include all covariates (floral abundance, floral richness, BB abundance, BB richness, transect number) in the models, understanding that we have collinearity problems, and see which ones explain the most. Because of the (curvilinear) collinearity between elevation and the other covariates, though, this approach would make it harder to detect and interpret the infleunce of elevation, which is the main question we’re interested in.

  2. The second approach is to say that while network metrics may, indeed, be influenced by abundance and richness covariates, these are just mediated effects of elevation. We may not care exactly how (i.e. via which mediating covariate) elevation influences network metrics; we might just want to know the shape of the response. But I think we must include transects as a covariate even in this approach; otherwise, we would need to go down the road of bootstrapping/rarefaction, and I think we all agree that we don’t want to do that.

My recommendation would be to go with option 2, and that is what will be presented below.

Network/group-level metrics

For these models, we will use the bs = “fs” factor smooth method. This preserves degrees of freedom by forcing all smooths to have the same wiggliness.

Results

  • Niche overlap (15% explained)
    • Linear negative effect of elevation
    • No significant effect of transect number
  • C score (30% explained)
    • Positive quadratic/cubic effect of elevation
    • No significant effect of transect number

Discussion These models have a lot of unexplained variation, but they suggest that the diets of co-existing bumble bee species become increasingly dissimilar (niche overlap) or disaggregated (C score) with increasing elevation. This pattern does not appear to be influenced significantly by transect count.

# Web asymmetry
gam_web.asym <- gam(web.asymmetry ~ 
                            s(elev.mean, year, bs = "fs", k = 5) +
                            s(bb.transects, year, bs = "fs", k = 4),
                    data = netmet,
                    family = "gaussian",
                    method = "REML",
                    select = FALSE) %>% getViz()
          
check.gamViz(gam_web.asym) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-1.609747e-05,4.434315e-05]
## (score -42.69759 & scale 0.01290583).
## Hessian positive definite, eigenvalue range [1.090476e-06,36.6072].
## Model rank =  28 / 28 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                            k'      edf k-index p-value    
## s(elev.mean,year)    1.50e+01 1.21e-04    0.74  <2e-16 ***
## s(bb.transects,year) 1.20e+01 5.52e+00    1.14    0.89    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(gam_web.asym)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## web.asymmetry ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects, 
##     year, bs = "fs", k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.50578    0.06483  -7.801 5.27e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df     F  p-value    
## s(elev.mean,year)    0.0001214     14  0.00    0.544    
## s(bb.transects,year) 5.5163473     10 10.58 7.83e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.592   Deviance explained = 62.3%
## -REML = -42.698  Scale est. = 0.012906  n = 74
plot(sm(gam_web.asym, 1)) +
  l_fitLine(alpha = 0.6)

plot(sm(gam_web.asym, 2)) +
  l_fitLine(alpha = 0.6)

# Niche overlap
gam_niche.overlap <- gam(niche.overlap.HL ~ 
                            s(elev.mean, year, bs = "fs", k = 5) +
                            s(bb.transects, year, bs = "fs", k = 4),
                         data = netmet,
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()
          
check.gamViz(gam_niche.overlap) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 13 iterations.
## Gradient range [-4.387283e-06,8.27703e-06]
## (score -52.89884 & scale 0.01218919).
## Hessian positive definite, eigenvalue range [2.061461e-06,36.53777].
## Model rank =  28 / 28 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                            k'      edf k-index p-value  
## s(elev.mean,year)    1.50e+01 2.32e+00    0.83   0.065 .
## s(bb.transects,year) 1.20e+01 2.56e-05    1.06   0.610  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(gam_niche.overlap)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## niche.overlap.HL ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects, 
##     year, bs = "fs", k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.28475    0.01284   22.18   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df    F p-value   
## s(elev.mean,year)    2.321e+00     14 0.73 0.00668 **
## s(bb.transects,year) 2.558e-05     11 0.00 0.86758   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.123   Deviance explained = 15.1%
## -REML = -52.899  Scale est. = 0.012189  n = 74
plot(sm(gam_niche.overlap, 1)) +
  l_fitLine(alpha = 0.6)

plot(sm(gam_niche.overlap, 2)) +
  l_fitLine(alpha = 0.6)

# C score
gam_C.score <- gam(C.score.HL ~ 
                            s(elev.mean, year, bs = "fs", k = 5) +
                            s(bb.transects, year, bs = "fs", k = 4),
                   data = netmet,
                   family = "gaussian",
                   method = "REML",
                   select = FALSE) %>% getViz()
          
check.gamViz(gam_C.score) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 12 iterations.
## Gradient range [-5.627583e-06,3.739096e-06]
## (score -36.37735 & scale 0.01741642).
## Hessian positive definite, eigenvalue range [2.064089e-06,36.63284].
## Model rank =  28 / 28 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                          k'    edf k-index p-value  
## s(elev.mean,year)    15.000  4.934    0.83    0.04 *
## s(bb.transects,year) 12.000  0.697    1.09    0.76  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(gam_C.score)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## C.score.HL ~ s(elev.mean, year, bs = "fs", k = 5) + s(bb.transects, 
##     year, bs = "fs", k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.43572    0.04316    10.1 4.09e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                         edf Ref.df     F  p-value    
## s(elev.mean,year)    4.9343     14 1.379 0.000469 ***
## s(bb.transects,year) 0.6967     10 0.134 0.075160 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.238   Deviance explained = 29.7%
## -REML = -36.377  Scale est. = 0.017416  n = 74
plot(sm(gam_C.score, 1)) +
  l_fitLine(alpha = 0.6)

plot(sm(gam_C.score, 2)) +
  l_fitLine(alpha = 0.6)

Species-level d’

First let’s look at the abundance and distribution of each species in each year. It’s probably not worth trying to fit the model to species that are extremely rare or limited in their distirbution. We will drop B. humilis and B. lapidarius from 2010, B. humilis, B. gerstaeckeri, and B. mucidus from 2011, and B. humilis, B. hypnorum, and B. mucidus from 2012.

bb_range.year <- net %>%
  left_join(site_data) %>%
  group_by(site, date, elev.mean, bb.sp, year) %>%
  summarize(abund = n()) %>%
  group_by(site, elev.mean, bb.sp, year) %>%
  summarize(abund = mean(abund)) %>%
  group_by(bb.sp) %>%
  mutate(elev.floor = min(elev.mean),
         elev.ceiling = max(elev.mean),
         elev.range = elev.ceiling - elev.floor,
         elev.med = elev.floor + ((elev.ceiling - elev.floor)/2))

ggplot(bb_range.year, aes(reorder(bb.sp, elev.med), elev.mean)) +
  geom_line(size = 2, alpha = 0.25) +
  geom_point(aes(size = abund), alpha = 0.5) +
  facet_wrap(~year, ncol = 1) +
  theme_light(12) +
  scale_size_continuous(name = "Abundance") +
  scale_color_discrete(name = "Abund. class") +
  xlab("BB species") +
  ylab("Elevation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

For species-level d’, we will use the “by = factor” factor smooth method. This allows different wiggliness and independent significance tests for each species, which I think is the right way to handle species-level differences in the response of d’ to elevation. We will also create separate models for each year so that we can drop the rare species from each year and avoid a complicated interaction term (bb.sp * year) in the by=factor call.

Results (only significant terms shown)

  • 2010 (33% explained)
    • B. gerstaeckeri (positive, linear)
    • B. pratorum (positive, quadratic)
    • B. pyrenaeus (negative, quadratic)
    • B. soroensis (positive, quadratic)
  • 2011 (41% explained)
    • B. hortorum (positive, linear)
    • B. psithyrus (positive, linear)
    • Transects (positive, quadratic)
  • 2012 (37% explained)
    • B. sensu-strictu (positive, cubic)
    • B. hortorum (negative, linear)
    • B. mendax (negative, leveling off above 1500)
    • B. pascuorum (negative, linear)
    • B. pyrenaeus (negative, linear)

Discussion The majority of species do not exhibit a significant relationship between elevation and discrimination, and which species do varies markedly by year. Only B. hortorum and B. pyrenaeus show significant relationships in more than a single year. Transect number was significant only in 2011. The most interesting finding is that the relationship between elevation and discrimination is generally positive in 2010 and 2011, but it changes to negative in 2012.

### Species-specific models
gam_d.2010 <- gam(d ~ 
                    bb.sp +
                    s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                    s(bb.transects, k = 4),
                  data = filter(spmet, year == 2010 & 
                                  !bb.sp %in% c("humi", "lapi")),
                  family = "gaussian",
                  method = "REML",
                  select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2010) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 16 iterations.
## Gradient range [-1.237157e-06,5.174296e-06]
## (score 31.42053 & scale 0.05353188).
## Hessian positive definite, eigenvalue range [4.464382e-07,83.51918].
## Model rank =  73 / 73 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                          k'  edf k-index p-value  
## s(elev.mean):bb.spbss  4.00 1.00    0.92   0.140  
## s(elev.mean):bb.spgers 4.00 1.00    0.92   0.145  
## s(elev.mean):bb.sphort 4.00 1.00    0.92   0.120  
## s(elev.mean):bb.sphypn 4.00 1.00    0.92   0.145  
## s(elev.mean):bb.spjone 4.00 1.00    0.92   0.085 .
## s(elev.mean):bb.spmend 4.00 1.00    0.92   0.115  
## s(elev.mean):bb.spmont 4.00 1.00    0.92   0.095 .
## s(elev.mean):bb.spmuci 4.00 1.36    0.92   0.150  
## s(elev.mean):bb.sppasc 4.00 1.00    0.92   0.130  
## s(elev.mean):bb.spprat 4.00 2.63    0.92   0.125  
## s(elev.mean):bb.sppsit 4.00 1.00    0.92   0.135  
## s(elev.mean):bb.sppyre 4.00 2.33    0.92   0.090 .
## s(elev.mean):bb.spsoro 4.00 2.35    0.92   0.115  
## s(elev.mean):bb.spwurf 4.00 1.00    0.92   0.100 .
## s(bb.transects)        3.00 1.00    0.90   0.080 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2010)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects, 
##     k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.35200    0.04632   7.599 2.25e-12 ***
## bb.spgers    0.05871    0.09037   0.650   0.5168    
## bb.sphort   -0.12442    0.07015  -1.774   0.0780 .  
## bb.sphypn   -0.01320    0.09957  -0.133   0.8947    
## bb.spjone   -0.01235    0.10967  -0.113   0.9105    
## bb.spmend   -0.24637    0.12576  -1.959   0.0518 .  
## bb.spmont    0.08340    0.15245   0.547   0.5851    
## bb.spmuci    0.08644    0.45429   0.190   0.8493    
## bb.sppasc   -0.10226    0.06966  -1.468   0.1440    
## bb.spprat    0.05830    0.06713   0.868   0.3865    
## bb.sppsit    0.15755    0.08034   1.961   0.0516 .  
## bb.sppyre    0.01964    0.15497   0.127   0.8993    
## bb.spsoro    0.09216    0.06646   1.387   0.1674    
## bb.spwurf    0.08295    0.06551   1.266   0.2072    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df     F p-value   
## s(elev.mean):bb.spbss  1.000  1.000 1.297 0.25639   
## s(elev.mean):bb.spgers 1.000  1.000 8.288 0.00452 **
## s(elev.mean):bb.sphort 1.000  1.000 1.784 0.18352   
## s(elev.mean):bb.sphypn 1.000  1.000 0.697 0.40518   
## s(elev.mean):bb.spjone 1.000  1.000 0.087 0.76805   
## s(elev.mean):bb.spmend 1.000  1.000 3.378 0.06789 . 
## s(elev.mean):bb.spmont 1.000  1.000 0.007 0.93355   
## s(elev.mean):bb.spmuci 1.358  1.593 0.065 0.84897   
## s(elev.mean):bb.sppasc 1.000  1.000 0.230 0.63212   
## s(elev.mean):bb.spprat 2.632  3.157 3.208 0.02442 * 
## s(elev.mean):bb.sppsit 1.000  1.000 2.882 0.09144 . 
## s(elev.mean):bb.sppyre 2.326  2.693 3.655 0.02507 * 
## s(elev.mean):bb.spsoro 2.345  2.849 2.965 0.04087 * 
## s(elev.mean):bb.spwurf 1.000  1.000 0.808 0.36996   
## s(bb.transects)        1.000  1.000 0.346 0.55741   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.196   Deviance explained =   33%
## -REML = 31.421  Scale est. = 0.053532  n = 196
print(plot(gam_d.2010, select = c(2, 10, 12, 13)), pages = 1)

gam_d.2011 <- gam(d ~ 
                     bb.sp +
                     s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                     s(bb.transects, k = 4),
                   data = filter(spmet, year == 2011 &
                                   !bb.sp %in% c("humi", "gers", "muci")),
                   family = "gaussian",
                   method = "REML",
                   select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2011) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-3.815379e-06,1.975298e-05]
## (score -22.71315 & scale 0.03124357).
## Hessian positive definite, eigenvalue range [5.346636e-07,104.0106].
## Model rank =  68 / 68 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                          k'  edf k-index p-value    
## s(elev.mean):bb.spbss  4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.sphort 4.00 1.24    0.64  <2e-16 ***
## s(elev.mean):bb.sphypn 4.00 1.34    0.64  <2e-16 ***
## s(elev.mean):bb.spjone 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.splapi 4.00 1.47    0.64  <2e-16 ***
## s(elev.mean):bb.spmend 4.00 1.93    0.64  <2e-16 ***
## s(elev.mean):bb.spmont 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.sppasc 4.00 1.24    0.64  <2e-16 ***
## s(elev.mean):bb.spprat 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.sppsit 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.sppyre 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.spsoro 4.00 1.00    0.64  <2e-16 ***
## s(elev.mean):bb.spwurf 4.00 2.02    0.64  <2e-16 ***
## s(bb.transects)        3.00 2.44    0.69  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2011)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects, 
##     k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.42663    0.03635  11.736  < 2e-16 ***
## bb.sphort   -0.09794    0.05313  -1.844  0.06670 .  
## bb.sphypn   -0.25762    0.06042  -4.264 3.08e-05 ***
## bb.spjone   -0.16543    0.05307  -3.117  0.00209 ** 
## bb.splapi   -0.16576    0.06144  -2.698  0.00756 ** 
## bb.spmend   -0.06793    0.09951  -0.683  0.49560    
## bb.spmont   -0.01154    0.05881  -0.196  0.84457    
## bb.sppasc   -0.02033    0.05378  -0.378  0.70575    
## bb.spprat   -0.08001    0.05182  -1.544  0.12417    
## bb.sppsit   -0.04304    0.05535  -0.777  0.43777    
## bb.sppyre   -0.01830    0.08372  -0.219  0.82719    
## bb.spsoro    0.06489    0.05140   1.262  0.20829    
## bb.spwurf    0.13937    0.05243   2.658  0.00847 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df      F  p-value    
## s(elev.mean):bb.spbss  1.000  1.000  0.104  0.74756    
## s(elev.mean):bb.sphort 1.241  1.441  7.304  0.00694 ** 
## s(elev.mean):bb.sphypn 1.339  1.599  1.047  0.43774    
## s(elev.mean):bb.spjone 1.000  1.000  1.831  0.17750    
## s(elev.mean):bb.splapi 1.472  1.789  1.856  0.24195    
## s(elev.mean):bb.spmend 1.926  2.212  2.326  0.10547    
## s(elev.mean):bb.spmont 1.000  1.000  1.828  0.17790    
## s(elev.mean):bb.sppasc 1.236  1.429  3.808  0.05709 .  
## s(elev.mean):bb.spprat 1.000  1.000  1.333  0.24958    
## s(elev.mean):bb.sppsit 1.000  1.000 19.677 1.48e-05 ***
## s(elev.mean):bb.sppyre 1.000  1.000  2.695  0.10221    
## s(elev.mean):bb.spsoro 1.000  1.000  2.279  0.13264    
## s(elev.mean):bb.spwurf 2.021  2.490  1.427  0.21996    
## s(bb.transects)        2.436  2.799  3.593  0.03613 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.319   Deviance explained = 40.8%
## -REML = -22.713  Scale est. = 0.031244  n = 235
print(plot(gam_d.2011, select = c(2, 10, 14)), pages = 1)

gam_d.2012 <- gam(d ~ 
                     bb.sp +
                     s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                     s(bb.transects, k = 4),
                   data = filter(spmet, year == 2012 &
                                   !bb.sp %in% c("humi", "hypn", "muci")),
                   family = "gaussian",
                   method = "REML",
                   select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2012) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 16 iterations.
## Gradient range [-1.188507e-06,8.313178e-06]
## (score -11.72191 & scale 0.03378739).
## Hessian positive definite, eigenvalue range [1.858343e-07,95.01981].
## Model rank =  68 / 68 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                          k'  edf k-index p-value  
## s(elev.mean):bb.spbss  4.00 3.10    0.89   0.020 *
## s(elev.mean):bb.spgers 4.00 1.00    0.89   0.070 .
## s(elev.mean):bb.sphort 4.00 1.13    0.89   0.060 .
## s(elev.mean):bb.spjone 4.00 1.00    0.89   0.050 *
## s(elev.mean):bb.splapi 4.00 1.61    0.89   0.045 *
## s(elev.mean):bb.spmend 4.00 1.83    0.89   0.030 *
## s(elev.mean):bb.spmont 4.00 1.00    0.89   0.055 .
## s(elev.mean):bb.sppasc 4.00 1.00    0.89   0.035 *
## s(elev.mean):bb.spprat 4.00 1.54    0.89   0.035 *
## s(elev.mean):bb.sppsit 4.00 1.00    0.89   0.050 *
## s(elev.mean):bb.sppyre 4.00 1.00    0.89   0.020 *
## s(elev.mean):bb.spsoro 4.00 2.31    0.89   0.060 .
## s(elev.mean):bb.spwurf 4.00 1.00    0.89   0.055 .
## s(bb.transects)        3.00 1.00    0.86   0.020 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2012)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ bb.sp + s(elev.mean, by = bb.sp, bs = "tp", k = 5) + s(bb.transects, 
##     k = 4)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.420864   0.036953  11.389  < 2e-16 ***
## bb.spgers    0.221879   0.096305   2.304  0.02234 *  
## bb.sphort   -0.085015   0.054129  -1.571  0.11799    
## bb.spjone   -0.211027   0.080103  -2.634  0.00914 ** 
## bb.splapi   -0.159219   0.076893  -2.071  0.03978 *  
## bb.spmend    0.210935   0.098327   2.145  0.03324 *  
## bb.spmont   -0.110492   0.070557  -1.566  0.11907    
## bb.sppasc   -0.024553   0.052687  -0.466  0.64175    
## bb.spprat   -0.077977   0.053391  -1.460  0.14585    
## bb.sppsit    0.115933   0.056859   2.039  0.04288 *  
## bb.sppyre    0.032636   0.088826   0.367  0.71373    
## bb.spsoro    0.016465   0.052236   0.315  0.75296    
## bb.spwurf   -0.003091   0.052170  -0.059  0.95282    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df     F p-value   
## s(elev.mean):bb.spbss  3.098  3.591 4.855 0.00272 **
## s(elev.mean):bb.spgers 1.000  1.000 0.355 0.55181   
## s(elev.mean):bb.sphort 1.134  1.254 5.310 0.01242 * 
## s(elev.mean):bb.spjone 1.000  1.000 0.188 0.66530   
## s(elev.mean):bb.splapi 1.614  1.961 0.344 0.68781   
## s(elev.mean):bb.spmend 1.830  2.198 3.308 0.03324 * 
## s(elev.mean):bb.spmont 1.000  1.000 0.491 0.48435   
## s(elev.mean):bb.sppasc 1.000  1.000 4.411 0.03704 * 
## s(elev.mean):bb.spprat 1.545  1.895 0.855 0.49405   
## s(elev.mean):bb.sppsit 1.000  1.000 0.694 0.40587   
## s(elev.mean):bb.sppyre 1.000  1.000 5.170 0.02411 * 
## s(elev.mean):bb.spsoro 2.308  2.807 1.306 0.39447   
## s(elev.mean):bb.spwurf 1.000  1.000 0.030 0.86283   
## s(bb.transects)        1.000  1.000 3.695 0.05608 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.263   Deviance explained = 37.1%
## -REML = -11.722  Scale est. = 0.033787  n = 217
print(plot(gam_d.2012, select = c(1, 3, 6, 8, 11)), pages = 1)

Let’s revisit option 1.

If we’re not committed to a story about elevation, can we explain our network metrics by including a larger suite of covariates, some of which may be collinear with elevation?

First, let’s inspect collinearity amongst the covariates. A PCA won’t be much help here since some of these relationships are nonlinear. So, we’ll make a pairs plot and interpret it manually. There’s a lot that could be said about these graphs, but here are some highlights.

For a start, we will use all our covariates except for floral transects, since that one has to most tangential relationship to d’ (it’s only effect would be via the estimation of floral species’ relative abundances). We will use the select = TRUE option to perform double-penalty variable selection; this potentially penalizes some terms out of the model.

Network/group-level metrics

For these models, we will use the bs = “fs” factor smooth method. This preserves degrees of freedom by forcing all smooths to have the same wiggliness.

Results

  • Niche overlap (37% explained)
    • Weakly quadratic negative effect of floral richness
    • Linear positive effect of (log) BB abundance
    • No significant effect of elevation
    • Linear or weakly quadratic positive effect of transect count in 2011 and 2012, but negative (weakly quadratic) effect in 2010.
  • C score (57% explained)
    • Cubic negative effect of BB richness.
    • Positive linear effect of elevation
    • Negative linear effect of transect number

Discussion These models explain a lot more variance than the elevation-only models did. Interestingly, the effect of elevation drops out entirely in the niche overlap model with the addition of the other covariates. This suggests that the effects of elevation may, indeed, be mediated by effects on abundance and diversity, but moreover that abundance and diversity affect network metrics apart from the influence of elevation.

The story for niche overlap seems to be that diets become more similar as floral richness goes down and bumble bee abundance goes up. The relationship with floral richness seems easy to interpret: when there are fewer options, bumble bees are forced to overlap more in their choices. The positive relationship with BB abundance might be a samping effect; when a lot of bumble bees are seen, it is more likely that their observed diets will overlap.

C score increased lineary with elevation, indicating that, as one moves upslope, BB diets tend to become more disaggregated (in a presence/absence sense – C score is basically an inverse and binary version of niche overlap). It also decreased linearly with transect number, which is probably a sampling effect: with higher sampling intensity, more links are observed, resulting in more observed overlap and less disaggregation. The negative cubic relationship with BB richness – and remember that BB richness and BB abundance are almost perfectly correlated – also suggests a sampling effect: when there are more and more kinds of BB, observed overlap is higher.

# Niche overlap
gam_niche.overlap.global <- gam(niche.overlap.HL ~ 
                                  s(bb.rich, k = 4) +
                                  s(fl.rich, k = 4) +
                                  s(log.fl.abund, k = 4) +
                                  s(log.bb.abund, k= 4) +
                                  s(elev.mean, year, bs = "fs", k = 5) +
                                  s(bb.transects, year, bs = "fs", k = 4),
                                data = netmet,
                                family = "gaussian",
                                method = "REML",
                                select = TRUE) %>% getViz()
          
check.gamViz(gam_niche.overlap.global) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-1.034418e-05,3.131483e-05]
## (score -55.64349 & scale 0.009635063).
## Hessian positive definite, eigenvalue range [3.24396e-07,36.04731].
## Model rank =  40 / 40 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                            k'      edf k-index p-value
## s(bb.rich)           3.00e+00 1.27e-05    0.97    0.32
## s(fl.rich)           3.00e+00 2.11e+00    1.11    0.78
## s(log.fl.abund)      3.00e+00 6.61e-01    1.05    0.55
## s(log.bb.abund)      3.00e+00 8.84e-01    0.89    0.13
## s(elev.mean,year)    1.50e+01 1.29e+00    0.94    0.28
## s(bb.transects,year) 1.20e+01 1.67e+00    1.05    0.67

summary(gam_niche.overlap.global)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## niche.overlap.HL ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, 
##     k = 4) + s(log.bb.abund, k = 4) + s(elev.mean, year, bs = "fs", 
##     k = 5) + s(bb.transects, year, bs = "fs", k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27462    0.01889   14.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df     F p-value    
## s(bb.rich)           1.268e-05      3 0.000 0.42032    
## s(fl.rich)           2.113e+00      3 4.819 0.00059 ***
## s(log.fl.abund)      6.610e-01      3 0.650 0.06306 .  
## s(log.bb.abund)      8.844e-01      3 2.550 0.00137 ** 
## s(elev.mean,year)    1.295e+00     14 0.158 0.13894    
## s(bb.transects,year) 1.671e+00     10 0.456 0.03922 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.311   Deviance explained = 37.4%
## -REML = -55.643  Scale est. = 0.0096351  n = 73
print(plot(gam_niche.overlap.global, allTerms = T), pages = 1)

plot(sm(gam_niche.overlap.global, 6)) +
  l_fitLine(alpha = 0.6) 

# C score
gam_cscore.global <- gam(C.score.HL ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           s(elev.mean, year, bs = "fs", k = 5) +
                           s(bb.transects, year, bs = "fs", k = 4),
                         data = netmet,
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()
          
check.gamViz(gam_cscore.global) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 15 iterations.
## Gradient range [-5.128725e-06,1.300864e-05]
## (score -37.88911 & scale 0.01162369).
## Hessian positive definite, eigenvalue range [2.399802e-06,34.07957].
## Model rank =  40 / 40 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                         k'   edf k-index p-value
## s(bb.rich)            3.00  2.78    1.20    0.94
## s(fl.rich)            3.00  2.34    1.02    0.52
## s(log.fl.abund)       3.00  1.60    0.95    0.33
## s(log.bb.abund)       3.00  1.14    1.13    0.87
## s(elev.mean,year)    15.00  1.93    1.08    0.68
## s(bb.transects,year) 12.00  1.26    1.09    0.73

summary(gam_cscore.global)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## C.score.HL ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, 
##     k = 4) + s(log.bb.abund, k = 4) + s(elev.mean, year, bs = "fs", 
##     k = 5) + s(bb.transects, year, bs = "fs", k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4809     0.0178   27.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df     F p-value  
## s(bb.rich)           2.781  2.953 4.604  0.0108 *
## s(fl.rich)           2.337  2.707 3.205  0.0680 .
## s(log.fl.abund)      1.603  1.962 2.295  0.1399  
## s(log.bb.abund)      1.138  1.246 0.321  0.5483  
## s(elev.mean,year)    1.935 14.000 0.362  0.0457 *
## s(bb.transects,year) 1.256 10.000 0.397  0.0313 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.49   Deviance explained = 56.8%
## -REML = -37.889  Scale est. = 0.011624  n = 73
print(plot(gam_cscore.global, allTerms = T), pages = 1)

plot(sm(gam_cscore.global, 4)) +
  l_fitLine(alpha = 0.6)

plot(sm(gam_cscore.global, 5)) +
  l_fitLine(alpha = 0.6) 

Species-level d’

Results (only significant terms shown)

Discussion This is interesting. I do not fully understand why, but the addition of abundance and diversity covariates actually strengthens the apparent effects of elevation. It could also have something to do with setting select = TRUE. The change is most evident in 2011, where most BB species now show a significant relationship with elevation. The overall patterns are the same, though. Elevation has a positive effect on discrimination in 2010 and 2011, but this switches to a negative relationship in 2012.

Conclusion

### Species-specific models
gam_d.2010.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           bb.sp +
                           s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2010 & 
                                         !bb.sp %in% c("humi", "lapi")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()

gam_d.2010.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           s(elev.mean, bb.sp, bs = "fs", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2010 & 
                                         !bb.sp %in% c("humi", "lapi")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2010.global) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-2.25218e-06,3.597569e-06]
## (score 20.06302 & scale 0.0569079).
## Hessian positive definite, eigenvalue range [9.160175e-07,91.72902].
## Model rank =  86 / 86 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'   edf k-index p-value
## s(bb.rich)          3.00  1.00    1.04    0.65
## s(fl.rich)          3.00  1.00    1.05    0.74
## s(log.fl.abund)     3.00  2.54    1.07    0.78
## s(log.bb.abund)     3.00  1.00    1.08    0.88
## s(elev.mean,bb.sp) 70.00  8.92    1.08    0.82
## s(bb.transects)     3.00  1.48    1.06    0.79
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2010.global)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) + 
##     s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) + 
##     s(bb.transects, k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.36559    0.02585   14.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df     F p-value   
## s(bb.rich)         1.000  1.000 0.019 0.89075   
## s(fl.rich)         1.000  1.000 0.000 0.99621   
## s(log.fl.abund)    2.541  2.852 1.542 0.16655   
## s(log.bb.abund)    1.000  1.000 0.522 0.47099   
## s(elev.mean,bb.sp) 8.918 66.000 0.285 0.00573 **
## s(bb.transects)    1.485  1.756 0.719 0.37097   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.118   Deviance explained = 19.3%
## -REML = 20.063  Scale est. = 0.056908  n = 189
print(plot(gam_d.2010.global, allTerms = TRUE), pages = 1)

print(plot(gam_d.2010.global, select = c(6, 7, 10, 14, 15, 16, 17)), pages = 1)

gam_d.2011.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           bb.sp +
                           s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2011 &
                                         !bb.sp %in% c("humi", "gers", "muci")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()

gam_d.2011.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           s(elev.mean, bb.sp, bs = "fs", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2011 &
                                         !bb.sp %in% c("humi", "gers", "muci")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2011.global) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-5.129729e-05,0.0001512141]
## (score -52.69027 & scale 0.02657087).
## Hessian positive definite, eigenvalue range [2.727704e-06,114.903].
## Model rank =  81 / 81 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'   edf k-index p-value    
## s(bb.rich)          3.00  1.00    0.79  <2e-16 ***
## s(fl.rich)          3.00  1.00    0.79  <2e-16 ***
## s(log.fl.abund)     3.00  2.20    0.79  <2e-16 ***
## s(log.bb.abund)     3.00  1.66    0.77  <2e-16 ***
## s(elev.mean,bb.sp) 65.00 18.68    0.79  <2e-16 ***
## s(bb.transects)     3.00  2.50    0.81  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2011.global)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) + 
##     s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) + 
##     s(bb.transects, k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.3782     0.0297   12.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df     F  p-value    
## s(bb.rich)          1.000  1.000 7.494  0.00672 ** 
## s(fl.rich)          1.000  1.000 1.204  0.27386    
## s(log.fl.abund)     2.204  2.612 2.502  0.12074    
## s(log.bb.abund)     1.663  2.014 1.541  0.22276    
## s(elev.mean,bb.sp) 18.682 64.000 1.589 3.91e-13 ***
## s(bb.transects)     2.501  2.810 5.521  0.00660 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.421   Deviance explained = 48.8%
## -REML = -52.69  Scale est. = 0.026571  n = 235
print(plot(gam_d.2011.global, allTerms = TRUE), pages = 1)

print(plot(gam_d.2011.global, select = c(1, 2, 5, 6, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18)), pages = 1)

gam_d.2012.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           bb.sp +
                           s(elev.mean, by = bb.sp, bs = "tp", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2012 &
                                         !bb.sp %in% c("humi", "hypn", "muci")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()

gam_d.2012.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           s(elev.mean, bb.sp, bs = "fs", k = 5) +
                           s(bb.transects, k = 4),
                         data = filter(spmet, year == 2012 &
                                         !bb.sp %in% c("humi", "hypn", "muci")),
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()

gam_d.2012.global <- gam(d ~ 
                           s(bb.rich, k = 4) +
                           s(fl.rich, k = 4) +
                           s(log.fl.abund, k = 4) +
                           s(log.bb.abund, k= 4) +
                           s(elev.mean, bb.sp, bs = "fs", k = 5) +
                           s(bb.transects, k = 4),
                         data = spmet,
                         family = "gaussian",
                         method = "REML",
                         select = FALSE) %>% getViz()
          
check.gamViz(gam_d.2012.global) 
## 
## Method: REML   Optimizer: outer newton
## full convergence after 17 iterations.
## Gradient range [-1.927926e-05,9.303341e-05]
## (score -45.08819 & scale 0.04610161).
## Hessian positive definite, eigenvalue range [1.798701e-06,325.6101].
## Model rank =  96 / 96 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'   edf k-index p-value    
## s(bb.rich)          3.00  1.00    0.89  <2e-16 ***
## s(fl.rich)          3.00  2.11    0.89  <2e-16 ***
## s(log.fl.abund)     3.00  1.35    0.89  <2e-16 ***
## s(log.bb.abund)     3.00  1.00    0.88   0.005 ** 
## s(elev.mean,bb.sp) 80.00 11.84    0.87  <2e-16 ***
## s(bb.transects)     3.00  1.00    0.88  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

summary(gam_d.2012.global)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## d ~ s(bb.rich, k = 4) + s(fl.rich, k = 4) + s(log.fl.abund, k = 4) + 
##     s(log.bb.abund, k = 4) + s(elev.mean, bb.sp, bs = "fs", k = 5) + 
##     s(bb.transects, k = 4)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.37188    0.02279   16.32   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df     F  p-value    
## s(bb.rich)          1.000  1.000 0.021   0.8844    
## s(fl.rich)          2.114  2.522 2.536   0.0724 .  
## s(log.fl.abund)     1.346  1.612 6.837   0.0121 *  
## s(log.bb.abund)     1.000  1.000 0.305   0.5808    
## s(elev.mean,bb.sp) 11.841 77.000 0.870 1.21e-10 ***
## s(bb.transects)     1.000  1.000 3.882   0.0492 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.123   Deviance explained = 14.7%
## -REML = -45.088  Scale est. = 0.046102  n = 657
print(plot(gam_d.2012.global, select = c(1, 3, 4, 6, 10, 12, 15)), pages = 1)

plot(gam_d.2012.global, select = c(5))

plot(gam_d.2011.global, select = c(5))

plot(gam_d.2010.global, select = c(5))

Conclusion

When niche overlap and C score are modeled only as functions of elevation and transect count, we can find interpretable, if weak, elevational patterns. When we include the other covariates, though, the influence of elevation is largely eclipsed by that of BB abundance/richness, and I am not confident that these relationships can be interpreted as anything more than sampling effects. This casts doubt on the legitimacy of the elevation patterns found when abundance and richness covariates are excluded; it’s hard to make the case that they are not just sampling effects, too.

The preliminary analysis below are not finished yet, so I won’t include the results at this stage.

Module analysis

First for the metaweb

Diet similarity by ordination and cluster analysis

Preliminary visualizations

Network metrics ~ elevation

Group-level metrics ~ elevation

Species-level d’ ~ elevation

Species-level d’ ~ bb.sp, year, and tongue length